function Y = spectral_embedding(X, ndim, knnr, sigma, wflag)
    if nargin < 2 || isempty(ndim)
        ndim = 10;
    end
    if nargin < 3 || isempty(knnr)
        knnr = 0.075;
    end
    if nargin < 4 || isempty(sigma)
        sigma = 1;
    end
    if nargin < 5 || isempty(wflag)
        wflag = 'simple';
    end
    
    %%% get general parameters %%%
    np = size(X, 1);
    nr = length(knnr);
    Y = cell(1, nr);
    
    for ii = 1: nr
        %%% compute weight matrix & laplacian %%%
        knn = max(20, round(knnr(ii) * np));
        W = compute_geo_laplacian(X, knn, wflag, sigma);
        
        %%% compute laplacian & eigen-decomp matrix %%%
        D = diag(sum(W, 2));
        Dh = diag(sum(W, 2) .^ -0.5);
        Dph = diag(sum(W, 2) .^ 0.5);
        Lt = D - W;
        LL = Dh * Lt * Dh;
        
        %%% compute eigen decomposition %%%
        [v, d] = eigs(LL, ndim + 1, 'sa');
        
        %%% get embedding %%%
        Y{ii} = Dph * v(:, 2: ndim + 1);
%         Y = Y ./ sqrt(sum(Y .^ 2, 2));

        %%% update X %%%
        X = Y{ii};
    end
end





